%% 2023/02/15
% Figure 3 (anyon pairing with U = 0), case (a-d)

clear
close all


%% Directories containing the data

theta_list = 0:0.5:1; % Statistical phases, in unit of pi
N = numel(theta_list);

J = 2*pi * 10.5; % Tunneling amplitude

common_root = 'mat_files\m2023_01_';
filename_list{1} = strcat(common_root, '13_4_2D2_qw_n2_6Er_modfreq3_1010_phase0.mat');
filename_list{2} = strcat(common_root, '09_1_2D2_qw_n2_6Er_modfreq3_1010_phase0p5.mat');
filename_list{3} = strcat(common_root, '05_1_2D2_qw_n2_6Er_modfreq3_1010_phase1p0.mat');


%% Other things that can be useful

% Matlab default colors for plotting
c1 = [0 0.4470 0.7410];
c2 = [0.8500 0.3250 0.0980];
c3 = [0.9290 0.6940 0.1250];


%% Load the data 

for i = 1:N
    filename = filename_list{i};
    data = load(filename);
    density{i} = data.bo_compressed;
    correlations{i} = data.corr_compressed;
    time{i} = J * 10^(-3) *  data.uniquetvals ;
end


%% Plot the conditional density profiles

for i = 1:N

    % Extract size grid
    [N_t, M] = size(density{i});
    center = [(M+1)/2, (M+1)/2];
    x_list = (1:M) - M/2;
    [X1, X2] = meshgrid(x_list, x_list);

    % Conditional correlations
    correlations_close = zeros([N_t, M, M]);
    correlations_far = zeros([N_t, M, M]);
    
    % Condition on distance between particles
    D_close = 3;
    index_close = abs(X1 - X2) < D_close;
    index_far = ~index_close;
    
    for j = 1:N_t
        correlations_aux = correlations{i}(:,:,j);
        correlations_close(j,:,:) = correlations_aux .* index_close;
        correlations_far(j,:,:) = correlations_aux .* index_far;
    end
    
    % Resum to get conditional density profiles
    density_close{i} = sum(correlations_close, 3);
    density_far{i} = sum(correlations_far, 3);


    %% Conditional density profiles: Close

    plot_figure = 1;
    if plot_figure
        figure
        imagesc([12, -11], time{i}, density_close{i})
        xlabel('Position j (sites)')
        ylabel('\tau')
        clim([0 1.1])
        colorbar
        xticks([-10, 0, 10])
        yticks([0, 1, 2, 3])
        title(['\theta = ' num2str(theta_list(i))])
    end


    %% Conditional density profiles: Far

    plot_figure = 1;
    if plot_figure
        figure
        imagesc([12, -11], time{i}, density_far{i})
        xlabel('Position j (sites)')
        ylabel('\tau')
        xticks([-10, 0, 10])
        yticks([0, 1, 2, 3])
        cb = colorbar();
        clim([0 1.1])
        title(['\theta = ' num2str(theta_list(i))])
    end

end
